Линейная алгебра. Практическая работа №11
Программные средства¶
Поверхности второго порядка¶
Общий вид
$F(x,y,z) = 0$
Параметрические
$(t_1, t_2)$ - внутренние координаты поверхности
$x = X(t_1, t_2); y = Y(t_1, t_2); z = Z(t_1, t_2)$
Как функция двух переменных
$z = f(x,z)$
Эллиптический тип¶
Сфера
$x^2 + y^2 + z^2 = R^2$
Параметрические
$x = R \cos t_1 \sin t_2; y = R \sin t_1 \sin t_2; z = R \cos t_2$
$0 \leq t_1 < 2\pi, 0 \leq t_2 \leq \pi$
Эллипсоид
$\frac{x^2}{p^2} + \frac{y^2}{q^2} + \frac{z^2}{r^2} = 1$
Мнимый эллипсоид
$\frac{x^2}{p^2} + \frac{y^2}{q^2} + \frac{z^2}{r^2} = -1$
Мнимый конус
$\frac{x^2}{p^2} + \frac{y^2}{q^2} + \frac{z^2}{r^2} = 0$
In [1]:
var('x,y,z,R,p,q,r,t1,t2')
eq1=[x^2+y^2+z^2-R^2==0,
[x==R*cos(t1)*sin(t2),
y==R*sin(t1)*sin(t2),
z==R*cos(t2)],
x^2/p^2+y^2/q^2+z^2/r^2==1]
sol1=[solve([eq1[0],R==4],R,z),
solve([eq1[1][0],eq1[1][1],eq1[1][2],R==3],R,x,y,z),
solve([eq1[2],p==1,q==2,r==3],p,q,r,z)]
for s in sol1: pretty_print(s)
In [2]:
var('x,y,z,t1,t2'); s=(-4,4)
colors=['#3636ff','#ff3636','#ff36ff']
f1=[x^2+y^2+z^2-4^2,(3*cos(t1)*sin(t2),
3*sin(t1)*sin(t2),3*cos(t2)),
x^2/1^2+y^2/2^2+z^2/3^2-1]
s1='<left><font color="%s">$(x,y,z) = %s$</font></left>'
s2='<left><font color="%s">$%s = 0$</font></left>'
pretty_print(html(s1%(colors[1],latex(f1[1]))))
[pretty_print(html(s2%(colors[i],latex(f1[i])))) for i in [0,2]]
(implicit_plot3d(f1[0],s,s,s,opacity=.3,color=colors[0])+\
parametric_plot3d(f1[1],(0,2*pi),(0,pi),opacity=.5,color=colors[1])+\
implicit_plot3d(f1[2],s,s,s,plot_points=20,
color=((z).function(x,y,z),colormaps.bwr_r)))\
.show(frame=False)
Выбор цветовых решений.
In [3]:
for i in [0..18]: print(sorted(colormaps)[8*i:8*(i+1)])
In [4]:
import numpy,sympy,pylab,warnings
sympy.init_printing(use_unicode=True)
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.lines import Line2D
warnings.simplefilter("ignore", category=RuntimeWarning)
x,y,z,R,p,q,r,t1,t2=\
sympy.symbols('x,y,z,R,p,q,r,t1,t2')
eq1=[x^2+y^2+z^2-R^2,
[x-R*cos(t1)*sin(t2),
y-R*sin(t1)*sin(t2),
z-R*cos(t2)],
x^2/p^2+y^2/q^2+z^2/r^2-1]
sol1=[sympy.solve([eq1[0],R-4],R,z),
sympy.solve([eq1[1][0],eq1[1][1],eq1[1][2],R-3],R,x,y,z),
sympy.solve([eq1[2],p-1,q-2,r-3],p,q,r,z)]
for s in sol1: display(s); print('\n')
In [6]:
angle=numpy.linspace(0,2*numpy.pi,64)
T1,T2=numpy.meshgrid(angle,angle)
X=numpy.cos(T1)*numpy.sin(T2)
Y=numpy.sin(T1)*numpy.sin(T2)
Z=numpy.cos(T2)
labels=['$x^2+y^2+z^2=4^2$',
'$x=3cos(t_1)*sin(t_2)$'+\
'\n$y=3sin(t_1)*sin(t_2)$'+\
'\n$z=3cos(t_2)$',
'$x^2/1^2+y^2/2^2+z^2/3^2=1$']
colors=['#3636ff','#ff3636','#ff36ff']
fig=pylab.figure(figsize=(9,9))
ax=fig.add_subplot(111,projection='3d')
ax.plot_surface(4*X,4*Y,4*Z,color=colors[0],alpha=.1)
ax.plot_surface(3*X,3*Y,3*Z,color=colors[1],alpha=.1)
ax.plot_surface(1*X,2*Y,3*Z,cmap=pylab.cm.bwr_r,
linewidth=1,antialiased=True,
rstride=1,cstride=1)
fake2Dlines=[Line2D([0],[0],linestyle="none",
c=c,marker='o') for c in colors]
ax.legend(fake2Dlines,labels,loc=9); pylab.show()
In [7]:
%%r
X1<-Y1<-seq(from=-3.2,to=3.2,by=.9)
sphere1<-function(R,X,Y){
func<-function(x,y){(-x^2-y^2+R^2)^.5}
Z<-outer(X,Y,func); Z[is.na(Z)]<-0; list(Z,-Z)}
T1<-seq(from=0,to=2*pi,length.out=12)
T2<-seq(from=0,to=pi,length.out=12)
sphere2<-function(R,t1,t2){
list(R*cos(t1)*sin(t2),R*sin(t1)*sin(t2),R*cos(t2))}
X3<-seq(from=-.9,to=.9,by=.22); Y3<-2*X3
ellipsoid1<-function(p,q,r,X,Y){
func<-function(x,y){r*(-x^2/p^2-y^2/q^2+1)^.5}
Z<-outer(X,Y,func); Z[is.na(Z)]<-0; list(Z,-Z)}
Z1<-sphere1(4,X1,Y1); XYZ2<-sphere2(3,T1,T2)
Z3<-ellipsoid1(1,2,3,X3,Y3)
dZ1=data.frame(Z1[1])
rownames(dZ1)<-X1; colnames(dZ1)<-Y1; dZ1
Гиперболический тип¶
Двуполостный гиперболоид
$\frac{x^2}{p^2} + \frac{y^2}{q^2} - \frac{z^2}{r^2} = -1$
Однополостный гиперболоид
$\frac{x^2}{p^2} + \frac{y^2}{q^2} - \frac{z^2}{r^2} = 1$
Конус
$\frac{x^2}{p^2} + \frac{y^2}{q^2} - \frac{z^2}{r^2} = 0$
In [8]:
var('x,y,z,p,q,r')
eq2=[x^2/p^2+y^2/q^2-z^2/r^2==-1,
x^2/p^2+y^2/q^2-z^2/r^2==1,
x^2/p^2+y^2/q^2-z^2/r^2==0]
sol2=[solve([eq2[0],p==1,q==2,r==3],p,q,r,z),
solve([eq2[1],p==3,q==4,r==5],p,q,r,z),
solve([eq2[2],p==2,q==3,r==4],p,q,r,z)]
for s in sol2: pretty_print(s)
In [9]:
var('x,y,z'); s=(-7,7); ops=[.5,.7,.3]
colors=['#ff36ff','#ff3636','#3636ff']
f2=[x^2/1^2+y^2/2^2-z^2/3^2+1,
x^2/3^2+y^2/4^2-z^2/5^2-1,
x^2/2^2+y^2/3^2-z^2/4^2]
st='<left><font color="%s">$%s = 0$</font></left>'
[pretty_print(html(st%(colors[i],latex(f2[i]))))
for i in [0..2]]
sum([implicit_plot3d(f2[i],s,s,s,opacity=ops[i],mesh=True,
color=colors[i],plot_points=20)
for i in [0..2]]).show(frame=False)
In [10]:
x,y,z,p,q,r=sympy.symbols('x,y,z,p,q,r')
eq2=[x^2/p^2+y^2/q^2-z^2/r^2+1,
x^2/p^2+y^2/q^2-z^2/r^2-1,
x^2/p^2+y^2/q^2-z^2/r^2]
sol2=[sympy.solve([eq2[0],p-1,q-2,r-3],p,q,r,z),
sympy.solve([eq2[1],p-3,q-4,r-5],p,q,r,z),
sympy.solve([eq2[2],p-2,q-3,r-4],p,q,r,z)]
for s in sol2: display(s); print('\n')
In [11]:
X=numpy.linspace(-7,7,128); X,Y=numpy.meshgrid(X,X)
def f1(x,y): return 3/2*numpy.sqrt(4*x**2+y**2+4)
def f2(x,y): return 5/12*numpy.sqrt(16*x**2+9*y**2-144)
def f3(x,y): return 2/3*numpy.sqrt(9*x**2+4*y**2)
Z=[f1(X,Y),f2(X,Y),f3(X,Y),-f1(X,Y),-f2(X,Y),-f3(X,Y)]
for z in Z: z[z<-7]=numpy.NaN; z[z>7]=numpy.NaN
fig=pylab.figure(figsize=(9,9))
ax=fig.add_subplot(111,projection='3d')
colors=2*['#ff36ff','#ff3636','#3636ff']
alphas=2*[.6,.4,.2]
[ax.plot_wireframe(X,Y,Z[i],color=colors[i],
alpha=alphas[i]) for i in range(6)]
labels=['$x^2/1^2+y^2/2^2-z^2/3^2=-1$',
'$x^2/3^2+y^2/4^2-z^2/5^2=1$',
'$x^2/2^2+y^2/3^2-z^2/4^2=0$']
fake2Dlines=[Line2D([0],[0],linestyle="none",c=c,
marker='o') for c in colors[:3]]
ax.legend(fake2Dlines,labels,loc=9); pylab.show()
In [12]:
%%r
X1<-seq(from=-0.9,to=0.9,by=0.22); Y1<-2*X1
hyperboloid1<-function(p,q,r,X,Y){
func<-function(x,y){r*(x^2/p^2+y^2/q^2+1)^0.5}
Z<-outer(X,Y,func); list(Z,-Z)}
X2<-seq(from=-5,to=5,by=1.1); Y2<-4*X2/3
hyperboloid2<-function(p,q,r,X,Y){
func<-function(x,y){r*(x^2/p^2+y^2/q^2-1)^0.5}
Z<-outer(X,Y,func); Z[is.na(Z)]<-0; list(Z,-Z)}
X3<-seq(from=-1.9,to=1.9,by=0.5); Y3<-3*X3/2
cone1<-function(p,q,r,X,Y){
func<-function(x,y){r*(x^2/p^2+y^2/q^2)^0.5}
Z<-outer(X,Y,func); list(Z,-Z)}
Z1<-hyperboloid1(1,2,3,X1,Y1)
Z2<-hyperboloid2(3,4,5,X2,Y2)
Z3<-cone1(2,3,4,X3,Y3)
dZ3=data.frame(Z3[1])
rownames(dZ3)<-X3; colnames(dZ3)<-Y3; dZ3
Параболический тип¶
Эллиптический параболоид
$\frac{x^2}{p^2} + \frac{y^2}{q^2} = 2 * z$
Гиперболический параболоид
$\frac{x^2}{p^2} - \frac{y^2}{q^2} = 2 * z$
Эллиптический цилиндр
$\frac{x^2}{p^2} + \frac{y^2}{q^2} = 1$
Мнимый эллиптический цилиндр
$\frac{x^2}{p^2} + \frac{y^2}{q^2} = -1$
Гиперболический цилиндр
$\frac{x^2}{p^2} - \frac{y^2}{q^2} = 1$
Параболический цилиндр
$y^2 = 2 * p * x$
In [13]:
var('x,y,z,p,q')
eq3=[x^2/p^2+y^2/q^2==2*z,
x^2/p^2-y^2/q^2==2*z,
x^2/p^2+y^2/q^2-1==0,
x^2/p^2-y^2/q^2-1==0,
y^2==2*p*x]
sol3=[solve([eq3[0],p==1,q==2],p,q,z),
solve([eq3[1],p==2,q==3],p,q,z),
solve([eq3[2],p==3,q==4],p,q,y,z),
solve([eq3[3],p==2,q==3],p,q,y,z),
solve([eq3[4],p==2],p,y,z)]
for s in sol3: pretty_print(s)
In [14]:
var('x,y,z'); s=(-8,8); colors=['#ff3636','#3636ff']
f31=[x^2/1^2+y^2/2^2-2*z,x^2/2^2-y^2/3^2-2*z]
st='<left><font color="%s">$%s = 0$</font></left>'
[pretty_print(html(st%(colors[i],latex(f31[i]))))
for i in [0,1]]
sum([implicit_plot3d(f31[i],s,s,s,opacity=0.7,mesh=True,
plot_points=20,color=colors[i])
for i in [0,1]]).show()
In [15]:
var('x,y,z'); s=(-9,9)
colors=['#36ff36','#36ffff','#ff36ff']
f32=[x^2/3^3+y^2/4^2-1,x^2/2^3-y^2/3^2-1,y^2-2*2*x]
st='<left><font color="%s">$%s = 0$</font></left>'
[pretty_print(html(st%(colors[i],latex(f32[i]))))
for i in [0..2]]
sum([implicit_plot3d(f32[i],s,s,s,opacity=0.7,mesh=True,
plot_points=20,color=colors[i])
for i in [0..2]]).show()
In [16]:
x,y,z,p,q=sympy.symbols('x,y,z,p,q')
eq3=[x^2/p^2+y^2/q^2-2*z,
x^2/p^2-y^2/q^2-2*z,
x^2/p^2+y^2/q^2-1,
x^2/p^2-y^2/q^2-1,
y^2-2*p*x]
sol3=[sympy.solve([eq3[0],p-1,q-2],p,q,z),
sympy.solve([eq3[1],p-2,q-3],p,q,z),
sympy.solve([eq3[2],p-3,q-4],p,q,y,z),
sympy.solve([eq3[3],p-2,q-3],p,q,y,z),
sympy.solve([eq3[4],p-2],p,y,z)]
for s in sol3: display(s); print('\n')
In [18]:
X=numpy.linspace(-8,8,128)
X,Y=numpy.meshgrid(X,X)
def f1(x,y): return 1/2*x**2+1/8*y**2
def f2(x,y): return 1/8*x**2-1/18*y**2
Z=[f1(X,Y),f2(X,Y)]
colors=['#ff3636','#3636ff']; alphas=[0.8,0.4]
for z in Z: z[z<-8]=numpy.NaN; z[z>8]=numpy.NaN
fig=pylab.figure(figsize=(9,9))
ax=fig.add_subplot(111,projection='3d')
labels=['$x^2/1^2+y^2/2^2=2*z$',
'$x^2/2^2-y^2/3^2=2*z$']
[ax.plot_wireframe(X,Y,Z[i],color=colors[i],
alpha=alphas[i]) for i in range(2)]
fake2Dlines=[Line2D([0],[0],linestyle="none",
c=c,marker='o') for c in colors]
ax.legend(fake2Dlines,labels,loc=9); pylab.show()
In [136]:
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(X,Y,Z1, color='#ff3636', alpha=0.7)
ax.plot_wireframe(X,Y,Z2, color='#3636ff', alpha=0.5)
fake2Dline1 = mpl.lines.Line2D([0],[0], linestyle="none", c='#ff3636', marker = 'o')
fake2Dline2 = mpl.lines.Line2D([1],[1], linestyle="none", c='#3636ff', marker = 'o')
ax.legend([fake2Dline1, fake2Dline2],
['$x^2/1^2+y^2/2^2=2*z$', '$x^2/2^2-y^2/3^2=2*z$'], loc=9)
ax.set_xlim(-7,7); ax.set_ylim(-7,7); ax.set_zlim(-7,7)
plt.show()
In [20]:
X=numpy.linspace(-9,9,128)
X,Z=numpy.meshgrid(X,X)
def f3(x): return 4/3*(-x**2+9)**0.5
def f4(x): return 3/2*(x**2-4)**0.5
def f5(x): return -2*(x)**0.5
Y=[f3(X),f4(X),f5(X),-f3(X),-f4(X),-f5(X)]
colors=2*['#36ff36','#36ffff','#ff36ff']
alphas=2*[.5,.7,.3]
for y in Y: y[y<-9]=numpy.NaN; y[y>9]=numpy.NaN
fig=pylab.figure(figsize=(9,9))
ax=fig.add_subplot(111,projection='3d')
labels=['$x^2/3^2+y^2/4^2=1$',
'$x^2/2^2-y^2/3^2=1$',
'$y^2=2*2*x$']
[ax.plot_wireframe(X,Y[i],Z,\
color=colors[i],alpha=alphas[i])
for i in range(6)]
fake2Dlines=\
[Line2D([0],[0],linestyle="none",c=c,marker='o')
for c in colors[:3]]
ax.legend(fake2Dlines,labels,loc=9); pylab.show()
In [21]:
%%r
X1<-Y1<-X2<-Y2<-seq(from=-9,to=9,by=1)
X3<-Z3<-X4<-Z4<-Y5<-Z5<-seq(from=-5,to=5,by=0.1)
paraboloid1<-function(p,q,X,Y){
outer(X,Y,function(x,y){(x^2/p^2+y^2/q^2)/2})}
paraboloid2<-function(p,q,X,Y){
outer(X,Y,function(x,y){(x^2/p^2-y^2/q^2)/2})}
Z1<-paraboloid1(1,2,X1,Y1)
Z2<-paraboloid2(2,3,X2,Y2)
In [22]:
%%r
cylinder1<-function(p,q,X,Z){
for (x in X){for (z in Z){
y<-q*(1-x^2/p^2)^0.5
v<-c(x,y,z,x,-y,z)
if (sum(is.na(v))==0){
write.table(matrix(v,nrow=2,ncol=3,byrow=T),
file ="cylinder1.csv",append=T,
quote=F,col.names=F,row.names=F)}}}}
cylinder2<-function(p,q,X,Z){
for (x in X){for (z in Z){
y<-q*(x^2/p^2-1)^0.5
v<-c(x,y,z,x,-y,z)
if (sum(is.na(v))==0){
write.table(matrix(v,nrow=2,ncol=3,byrow=T),
file="cylinder2.csv",append=T,
quote=F,col.names=F,row.names=F)}}}}
cylinder3<-function(p,Y,Z){
for (y in Y){for (z in Z){
x<-y^2/2/p; v<-c(x,y,z)
write.table(matrix(v,nrow=1,ncol=3,byrow=T),
file ="cylinder3.csv",append=T,
quote=F,col.names=F,row.names=F)}}}
cylinder1(3,4,X3,Z3)
cylinder2(2,3,X4,Z4)
cylinder3(2,Y5,Z5)
"files are created"
In [24]:
files=['cylinder1.csv','cylinder2.csv',
'cylinder3.csv']
colors=['#36ff36','#36ffff','#ff36ff']
labels=['$x^2/3^2+y^2/4^2=1$',
'$x^2/2^2-y^2/3^2=1$',
'$y^2=2*2*x$']
XYZ=[numpy.array(numpy.loadtxt(files[i],\
delimiter=' ',unpack=True)) for i in range(3)]
fig=pylab.figure(figsize=(9,9))
ax=fig.add_subplot(111,projection='3d')
[ax.scatter(XYZ[i][0],XYZ[i][1],XYZ[i][2],s=1,\
c=colors[i],marker='1') for i in range(3)]
fake2Dlines=[Line2D([0],[0],linestyle="none",\
c=c,marker='1') for c in colors]
ax.legend(fake2Dlines,labels,loc=9); pylab.show()
Плоскости¶
Пара мнимых пересекающихся плоскостей
$\frac{x^2}{p^2} + \frac{y^2}{q^2} = 0$
Пара пересекающихся плоскостей
$\frac{x^2}{p^2} - \frac{y^2}{q^2} = 0$
Пара параллельных плоскостей
$y^2 = p^2$
Пара мнимых параллельных плоскостей
$y^2 = -p^2$
Пара совпадающих плоскостей
$y^2 = 0$
In [25]:
var('x,y,z,p,q')
eq4=[[x^2/p^2-y^2/q^2==0,p==1,q==2],
[y^2==p^2,p==1],y^2==0]
for e in eq4[:2]: pretty_print(solve(e,p,q,y))
s=(-3,3); f4=[x^2/1^2-y^2/2^2,y^2-1^2]
colors=['#ff3636','#ff36ff']
st='<left><font color="%s">$%s = 0$</font></left>'
[pretty_print(html(st%(colors[i],latex(f4[i]))))
for i in [0,1]]
(implicit_plot3d(f4[0],(x,-3,3),(y,-3,3),(z,-3,3),
color=(abs(sin(x*y*z)),colormaps.bwr))+\
implicit_plot3d(f4[1],s,s,s,
color=colors[1],opacity=0.7)).show()
In [27]:
x,y,z,p,q = sympy.symbols('x,y,z,p,q')
eq4=[[x^2/p^2-y^2/q^2,p-1,q-2],[y^2-p^2,p-1],y^2]
for e in eq4: display(sympy.solve(e,y)); print('\n')
X=numpy.linspace(-3,3,128)
X,Z=numpy.meshgrid(X,X)
def f1(x,z): return 2*x
def f2(x,z): return 1
Y=[f1(X,Z),-f1(X,Z),f2(X,Z),-f2(X,Z)]
for y in Y[:2]: y[y<-3]=numpy.NaN; y[y>3]=numpy.NaN
colors=['#ff3636','#ff36ff']
labels=['$x^2/1^2=y^2/2^2$','$y^2=1^2$']
fig=pylab.figure(figsize=(9,9))
ax=fig.add_subplot(111,projection='3d')
[ax.plot_surface(X,y,Z,cmap=pylab.cm.bwr,
linewidth=1,antialiased=True,
rstride=3,cstride=3)
for y in Y[:2]]
[ax.plot_wireframe(X,y,Z,color=colors[1],alpha=0.3)
for y in Y[2:]]
fake2Dlines=[Line2D([0],[0],linestyle="none",
c=c,marker='o') for c in colors]
ax.legend(fake2Dlines,labels,loc=9); pylab.show()
In [28]:
%%r
X<-Z<-seq(from=-3,to=3,by=0.1)
planes1<-function(p,q,X,Z){
for (x in X){for (z in Z){
y<-x*q/p
if (min(X)<=y&&y<=max(X)){
v<-c(x,y,z,x,-y,z)
v<-replace(v,is.na(v),0)
write.table(matrix(v,nrow=2,ncol=3,byrow=T),
file="planes1.csv",append=T,quote=F,
col.names=F,row.names=F)}}}}
planes2<-function(p,X,Z){
for (x in X){for (z in Z){
v<-c(x,p,z,x,-p,z)
write.table(matrix(v,nrow=2,ncol=3,byrow=T),
file="planes2.csv",append=T,
quote=F,col.names=F,row.names=F)}}}
planes1(1,2,X,Z); planes2(1,X,Z); "files are created"
In [29]:
X1,Y1,Z1=\
numpy.loadtxt('planes1.csv',delimiter=' ',unpack=True)
X2,Y2,Z2=\
numpy.loadtxt('planes2.csv',delimiter=' ',unpack=True)
fig=pylab.figure(figsize=(9,9))
ax=fig.add_subplot(111,projection='3d')
ax.scatter(X1,Y1,Z1,s=2,c=Z1,
cmap=pylab.cm.bwr,marker='1')
ax.scatter(X2,Y2,Z2,s=2,c='#ff36ff',marker='1')
labels=['$x^2/1^2=y^2/2^2$','$y^2=1^2$']
colors=['#ff3636','#ff36ff']
fake2Dlines=[Line2D([0],[0],linestyle="none",
c=c,marker='1') for c in colors]
ax.legend(fake2Dlines,labels,loc=9); pylab.show()
Задание.¶
Инициировать три поверхности второго порядка с помощью функций и
выполнить их трехмерные изображения,
используя SageMath, Python, R.